On the critical behavior of the one dimensional diffusive pair contact process 
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The phase transition of the one-dimensional, diffusive pair contact process (PCPD) is investigated 
by N cluster mean-field approximations and high precision simulations. The iV = 3, 4 cluster 
approximations exhibit smooth transition line to absorbing state by varying the diffusion rate D 
with /?2 = 2 mean-field order parameter exponent of the pair density. This contradicts with former 
N = 2 results, where two different mean-field behavior was found along the transition line. Extensive 
dynamical simulations on L = 10 s lattices give estimates for the order parameter exponents of the 
particles for 0.05 < D < 0.7. These data can support former two distinct class findings. However 
the gap between low and high D exponents is narrower than estimated previously and the possibility 
for interpreting numerical data as a single class behavior with exponents a — 0.21(1), /3 = 0.41(1) 
assuming logarithmic corrections is shown. Finite size scaling and cluster simulation results are also 
presented. 
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I. INTRODUCTION 

The exploration of nonequilibrium universality classes 
is current interest of research. In this area most sys- 
tems investigated exhibit phase transitions to absorbing 
states with such weak fluctuations from which no return 
is possible [1,2]. For a long time only the robust directed 
percolation (DP) universality class has been known [3,4]. 
Later systems with extra conservation laws and symme- 
tries were shown to belong to other universality classes 
[5-8]. In the past few years it turned out that there are 
novel classes in low dimensional reaction-diffusion sys- 
tems where neither classical bosonic field theory nor sym- 
metry arguments can give better understanding of the 
critical behavior [9]. This is probably due to the fact 
that in low dimensions topological constraints become 
effective, blocking the motion of reacting particles [10]. 
While bosonic field theories can not capture this feature, 
fermionic field theories have not been successful for such 
systems so far. In fact the critical behavior of such mod- 
els split according to fermionic or bosonic particles are 
involved in [11-14]. 

Recently novel universal behavior is reported in some 
low-dimensional reaction-diffusion models featured by 
production at pairs and single particle diffusion [13-21]. 
In these systems the production compete with pair anni- 
hilation and diffusion. If production wins steady states 
with finite particle density appear in (fermionic) models 
with hard-core repulsion, while in unrestricted (bosonic) 
models the density diverges. By lowering the produc- 
tion/annihilation rate a doublet of absorbing states with- 
out symmetries emerges. One of such states is completely 
empty, the other possesses a single wandering particle. 
In case of fermionic systems the transition to absorbing 
states is continuous with novel, yet not completely settled 
critical behavior. 

The field theory [13] describing bosonic particles could 
not be solved by standard renormalization procedures, 



but hinted at a transition with non-DP behavior. At 
the transition point of the Id model it predicts a density 
decay of the form 



p(t,Pc) oc 



ln(t) 



1/2 



(1) 



while in the inactive phase: p(t,p c ) oc i -1 / 2 . These were 
confirmed by simulations [10]. In case of fermionic par- 
ticles of this model (PCPD) density matrix renormaliza- 
tion group analysis [14], coherent anomaly extrapolation 
[16] and simulations [15,16] found novel kind of critical 
phase transition. However the critical exponents seem 
to depend on the diffusion strength D and different in- 
terpretations of data have been born. These embrace 
the possibilities of continuously changing exponents, two- 
universality classes [16] and single class with huge correc- 
tions [14,22]. 

Very recently well denned set of critical exponents 
are reported in different versions of binary production 
PCPD-likc processes [23]. However these simulations 
were done at a fixed, high diffusion/annihilation rate and 
as will be shown in Sect. IV the exponent estimates 
agree well with those of this paper in the high diffusion 
region. Even more recently two studies [24,25] reported 
non-universality in the dynamical behavior of the PCPD. 
While the former one by Dickman and Menezes explored 
different sectors (a reactive and a diffusive one) in the 
time evolution and gave nontrivial exponent estimates, 
the latter one by Hinrichsen provided a hypothesis that 
the ultimate long time behavior should be characterized 
by DP behavior. 

Just before the submission of this paper a preprint by 
Kockelkoren and Chate [26] showed extensive simulation 
results for a modified version of PCPD that is in be- 
tween fermionic and bosonic models. That means that 
they discard the single particle occupation constraint on 
the lattice but suppress multiple occupancy by an ex- 
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ponentially decreasing creation probability (p N ^ 2 ) of the 
particle number. They claim that their stochastic cellu- 
lar automaton (SCA) model shows smaller corrections to 
scaling than the PCPD and exhibit a single universality 
class transition. 

The two universality class scenario was backed by pair 
mean-field approximation [14] that showed two differ- 
ent mean-field behavior by varying D and simulations 
[16] for the order parameter density exponents. Such 
kind of mean-field behavior is absent if we replace the 
annihilation AA — > with coagulation AA — ► A [18]. 
By the investigation of the parity conserving version of 
the PCPD the mean-field and pair-mean-field approxi- 
mations resulted in similar phase diagram, but higher 
order cluster mean-field showed a single mean-field class 
behavior [21] and the authors concluded that the for ap- 
propriate description of such binary production models 
at least N — 3 clusters are needed. That mean-field be- 
havior was indeed found in d — d c — 2 by simulations 

In the present work I show N = 3, 4 cluster mean- 
field results for the PCPD model that again suggest sin- 
gle mean-field universality class. This does not necessar- 
ily imply that below d c = 2 only one class would exist. 
Higher precision simulations than that of [16] are also 
presented in the second part of this paper that provide 
better exponent estimates but still leave this question 
open. A single universality class scenario may be ac- 
cepted only if we assume logarithmic correction to data. 
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FIG. 1. Schematic phase diagram of the Id PCPD model. 
Circles correspond to simulation and DMRG results, solid line 
to site mean-field (N = 1), dashed line to pair-approximation 
(N = 2). Dot-dashed line shows N = 3, long-dashed N = 4 
cluster mean-field results discussed in Sect. III. 



According to pair mean-field approximations the phase 
diagram can be separated into two regions (see Fig.l). 
While for D > 1/7 the pair approximation gives the same 
p c (D) and exponents as the site mean-field, for D < 1/7-s 
the transition line breaks and the exponents are different 

a = 1, a 2 = 1, [3 = 1, fc = l. (6) 



II. THE PCPD MODEL 

A PCPD like binary spreading process was introduced 
in an early work by Grassberger [27] . Its preliminary sim- 
ulations in Id showed a non-DP type transition, but these 
results have been forgotten for a long time. The diffu- 
sive pair contact process (PCPD) introduced by Carlon 
et al [14] is controlled by two independent parameters: 
the probability of pair annihilation p and the probability 
of particle diffusion D. The dynamical rules are 



AA®, %AA AAA 
AA — > 00 
A0 «-» $A 



with rate(l -p)(l- D)/2 

with rate p(l — D) 

with rate D . (2) 



The site mean-field approximation gives a continuous 
transition at p — 1/3. For p < p c (D) the particle and 
pair densities exhibit singular behavior: 



p{oo,p) oc (p c -p) p p 2 (oo,p) tx (p c -pY 2 
while at p = p c {D) they decay as: 

p(t,p c ) <xt~ a , p 2 (t,p c ) «i~ Q2 , 
with the exponents: 

a = 1/2, a 2 = 1, 0=1, [3 2 = 2. 



(3) 
(4) 
(5) 



In the entire inactive phase the decay is characterized by 
the exponents: 



a 



1. 



Ct2 



(7) 



III. 



CLUSTER MEAN-FIELD RESULTS FOR 
PCPD 



The generalized, cluster mean-field approximation in- 
troduced by [28,29] was applied for the dynamical rules 
(2) of the Id fermionic lattice model. Master equations 
for N = 1, 2, 3, 4 block probabilities were set up 



dP N ({sj}) 
dt 



= f(PN({ Si })) 



(8) 



,A. The 



where site variables may take values: Sj ; 
equations could be solved numerically for the 9PjV ^ s 'D _ 
steady state condition. Taking into account spatial re- 
flection symmetries of Pjv({sj}) this involves 10 indepen- 
dent variables in case of N = 4. The particle {pip, D)) 
and pair {p 2 {p, D)) densities were expressed by Pv({si}) 
and the phase transition point p c (D) was located for sev- 
eral values of D. At p c (D) quadratic fitting of the form 
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a(p- Pc (D)) + b(p-p c (D)) 2 (9) 

was applied for p(p,D) and p2(p,D). The N = 1 and 2 
solutions reproduced the results of [14] for particle and 
pair densities. For N = 2 the two regions, correspond- 
ing to different leading order singularity of P2(p, D) with 
p2 = 1,2 were located by least square fit with the form 
(9). For N = 3,4 approximations smooth p c (d) phase 
transition lines are determined shown on Fig.l and tab- 
ulated in Table I. The quadratic fitting (9) resulted in 
leading order singularities (3=1 for particles and (32 — 2 
for pairs everywhere. These are in contradiction with 
the N — 2 approximation results similarly to the par- 
ity conserving binary process model case [21]. For that 
model simulations in 2d strengthened the single mean- 
field class behavior along p c (D) and it was conjectured 
that the pair approximation is an odd one. Here again 
I conclude that at least N > 2 level of approximation is 
necessary to obtain a correct mean-field behavior. 

The single mean-field class property does not neces- 
sarily mean that below d c a single class behavior should 
occur all along the p c (D) transition line. For example in 
a similar model that exhibits an additional global particle 
number conservation [8] such situation was found. There- 
fore I investigated by extensive simulations this question. 



The number of particles (N p ) is followed and the time 
is updated by 1/N p following a reaction (throughout the 
whole paper the time is measured by Monte Carlo steps 
(MCS)). The initial conditions were random distribution 
of particles with occupation probability 0.5. 

It was suggested in [24] that one may get smaller cor- 
rections to scaling if one excludes the purely diffusive sec- 
tor by averaging for states having at least one pair in the 
system. In the present simulations I did not find much 
effect (within statistical error margin) of such restrictions 
for the long time behavior. 5, 6. 




FIG. 3. Density decay times t 21 in Id PCPD at D = 0.5 
and p = 0.13351, 0.13352, 0.13353, 0.13356, 0.1336, 0.13363 
(top to bottom). 



IV. HIGH PRECISION SIMULATION RESULTS 

The simulations were performed on L = 10 5 sized rings 
with random sequential update version of PCPD evolv- 
ing by the following rules. A particle and a direction 
are selected randomly. One of the following reaction is 
performed: (a) a nearest neighbour exchange in the se- 
lected direction with probability D; (b) an annihilation 
with the nearest neighbour particle in the selected direc- 
tion with probability p(l — D); (c) a creation of a new 
particle in the selected direction at the second nearest 
neighbour empty site with probability (1 — p)(l — D) if 
the nearest neighbour is filled with a particle. 
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FIG. 2. Density decay times t 21 in Id PCPD at D = 0.7 
and p = 0.1574,0.15745,0.1575,0.15755,0.1576,0.1577 (top to 
bottom curves). 



A. Density decay simulations 

The critical point (p c ) for diffusion rates D — 0.05, 0.1, 
0.2, 0.5, 0.7 has been located by following the time evo- 
lution of the density decay. These simulations were done 
in two parts. First runs up to t max ~ 10 5 MCS and with 
high statistical averages (~ 10 4 ) were performed that al- 
lowed local slopes estimates of the density (p(t)) decay 
exponent a and p c . These simulations were extended by 
long time runs up to 10 7 - 10 8 MCS with 100 - 200 sam- 
ple numbers. The two sets of data are fitted together and 
are shown on Figs. 2, 3, 4, 
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FIG. 4. Density decay times t ' 21 in Id PCPD at D = 0.2 
and p = 0.111215,0.11217,0.11218,0.11219,0.1122 (top to bot- 
tom). 
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On all plots one can see up and down veering pit") 
curves in the long time limit - corresponding to active 
and absorbing phases - separated by a roughly straight 
line - corresponding to p c . As one can see for high dif- 
fusion rates (D > 0.2) scaling with exponent a ~ 0.21 
seems to a set for t >^ 3 x 10 4 MCS. This is in agree- 
ment with the first results provided for PCPD for high 
diffusion rates [16] and with the results of [17,24,26] for 
strong diffusions. 




t 



FIG. 5. Density decay times t 21 in Id PCPD at D = 0.1 
and p = 0.10686,0.10688,0.10689,0.1069,0.10691,0.10692 (top 
to bottom). 



In cases D = 0.05 and 0.1 straight lines on the log-log 
plot appear from t >~ 3 x 10 2 MCS with an exponent 
a = 0.245(5). This is in agreement with the results of 
[23] who considered the case when the coagulation and 
annihilation rate is three times the diffusion rate. This 
exponent is about 10% smaller than what was found in 
[16] but the two distinct class behavior seems to be sup- 
ported. 

0.95 , — — — — — — 



to scaling. For this reason I tried to apply logarithmic 
fitting to the data of the form 

p(t, Pc ) = [(a + b\n(t))/t] a . (10) 

One can find the corresponding exponents in Table II 
that are all in agreement with the value a = 0.21(1) in 
both the low and high diffusion regions. Here I applied 
least squares fitting for the most critical like curves such 
that the relative error in the sum of squares was at most 
0.0001. To confirm these results other critical exponents 
were investigated using the precise p c values of this sec- 
tion. 



B. Steady state simulations 

To estimate directly the order parameter exponent de- 
scribing the scaling 

p(oo, e) oc t p (11) 

off-critical, steady state densities had to be measured. 
Here again I used L = 10 5 system sizes. The density 
decay was followed for each D and ej = p c — pi values 
on logarithmic time scales until saturation effect was ob- 
served. Following that averaging of p{t) was done for 
about 100 samples within a time window that exceeds 
the saturation by a decade. I measured the effective ex- 
ponents like in [16] defined as 

ln(p(oo,ei)) - ln(p(oo,€i_i)) 



that are expected to converge to the true critical values 
in the e — > limit. 
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FIG. 6. Density decay times t - 21 in Id PCPD for D = 0.05 
and p = 0.10436, 0.10438,0.1044,0.10441,0.10442 (top to bot- 
tom). 



Although the upper critical dimension of PCPD is ex- 
pected to be at d c = 2 [21] one may not exclude the pos- 
sibility of a second critical dimension (d' c = 1) or topolog- 
ical effects in Id that may cause logarithmic corrections 
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FIG. 7. Effective j3 exponents for different diffusion rates. 
The circles correspond to D = 0.05, the squares to D = 0.1 
the diamonds to D — 0.2, the up-triangles to D = 0.5 and 
the down-triangles to D = 0.7. 



4 



As one can see on Fig. 7 the local slopes for D = 0.7 and 
D = 0.5 converge to (3 = 0.40(1) in agreement with the 
high diffusion rate results provided in [16]. This value 
is also close to Hinrichsen's estimate (0.38(6)) for the 
cyclically coupled model [17] and to Kockelkoren's value 
(0.37(2)) for the suppressed bosonic SCA model [26]. 

However for D = 0.05 and D — 0.1 extrapolations 
suggest [3 — 0.50(2). This is in agreement with Park's 
recent the results (~ 0.5) [23] but somewhat off the low- 
diffusion data of ref. [16] (0.57(2)) and from Dickman's 
estimates (0.55 — 0.45) [24]. The reason for these devi- 
ations is likely to be related to strong finite size effects, 
the complex way of scaling and the uncertainties of the 
p c values used. 



10' 
10 5 

J 

3 

e>T io 3 
J 

H 10' 
10"' 




100 1000 

L 



FIG. 8. Finite size scaling of tl (upper points) and pl(oo). 
The circles correspond to D = 0.05, the squares to D = 0.1 
the diamonds to D = 0.2, the up-triangles to D = 0.5 and the 
down-triangles to D = 0.7. The lines show power-law fittings 
applied for D = 0.7 data points. 



In case of the D — 0.2 curve one may observe an ex- 
trapolation to some intermediate value, but the curvature 
of the last points may also suggest a tendency towards 
the high D class data. Note that in earlier, lower scale 
simulations [16] the data for D = 0.2 showed low-D crit- 
ical behavior, strengthening the idea that some kind of 
very slow crossover happens here (although those results 
were obtained for a SCA version of PCPD). 

Similarly to dynamical simulations I tried the possibil- 
ity if logarithmic correction to scaling 

p(oo,e) = [e/(a + bln(e))f (13) 

could cure these "uncertainties" . As one can see in Ta- 
ble II the exponents for all D values satisfy scaling with 
(3 = 0.40(1) with logarithmic corrections of the form (13). 

C. Finite size scaling 

Finite size scaling at the critical point was performed 
for system sizes L = 32, 64, 128. ..1024. The quasi-steady 



state density (averaged over surviving samples) is ex- 
pected to scale according to 

p,(oo,p e ,L)«;L-M u± , (14) 

while the characteristic lifetime for half of the samples 
to reach the absorbing state scales with the dynamical 
exponent Z as 

t{ Pc ,L)kL z . (15) 

Since the system sizes are much smaller than in Sects. 
IV A and IV B one may expect stronger corrections to 
scaling. Indeed the power-law fitting for fi/v± results in 
values in the range 0.385 — 0.535 and for Z in the range 
1.75 — 2 depending on D. These results are shown of 
Fig. 8. Again the low- D data are in agreement with those 
of [14], [23] and [24], while the high-D data with those of 
[14], [26] and [17]. Just considering these ranges one can 
not distinguish this transition from the PC class (with 
/3/ux = 0.500(5) and Z = 1.75(1) [7]) that caused ini- 
tial debates in the literature [14-16]. Assuming single 
universality class corresponding to high- D data we may 
expect: 0/v x = 0.38(1) and Z = 1.75(15). 



V. CONCLUSIONS 

In this paper I addressed the long standing question of 
diffusion dependence of the phase transition of the PCPD 
model. The N = 3, 4 level cluster mean-field calculations 
confirmed a single mean-field universality class scenario 
similarly to the parity conserving version of this model 
[21]. Again the best conclusion one can draw from these 
data is that the N — 2 pair approximation is an odd one 
and we need at least N > 2 level of mean-filed to get the 
correct scaling behavior for binary production models. 

The extensive simulations have confirmed at least one 
set of the exponents - those for high diffusion- of the 
early results given in [16]. Data in the low diffusion 
range are in good agreement with other recent simula- 
tion results suggesting a different universality class. Al- 
though the scaling seems to set in much earlier in the 
low diffusion region than in the high diffusion range, a 
slow crossover to high-Z? behavior can be verified nu- 
merically assuming logarithmic corrections. Similar con- 
clusions can also be drawn from steady state simulation 
results. Although the two universality class picture pro- 
posed in [16] can not be excluded, data with logarithmic 
corrections may support a single class transition. Field 
theoretical arguments supporting or excluding logarith- 
mic corrections would be necessary. Note that in Id cou- 
pled systems logarithmic corrections are not rare at all. 

The finite size simulations could not give decisive sup- 
port for any of the the possible dependence on diffusion of 
this transition, but the range of results are in agreement 
with those of other numerical results of the literature. 
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Mean-filed exponents, the upper critical dimension and 
the lack of time reversal symmetry in this model seem 
to exclude the possibility of further crossovers to an 
ultimate DP critical behavior. Finally the insensitiv- 
ity to parity conservation in binary production models 
brings up the question of insensitivity for other conser- 
vation laws, hence binary production, diffusive models 
with global conservation might belong to the same class. 
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TABLE I. Summary of N = 2, 3, 4 approximation results 



D 


Pc 


P 


a 


0.7 


0.15745(1) 


0.39(1) 


0.214(5) 


0.5 


0.13353(1) 


0.414(16) 


0.206(7) 


0.2 


0.11218(1) 


0.402(8) 


0.217(8) 


0.1 


0.10688(1) 


0.407(7) 


0.206(7) 


0.05 


0.10439(1) 


0.411(10) 


0.216(9) 



TABLE II. Summary of simulation results assuming loga- 
rithmic corrections 
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